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Abstract 

An elementary Ising spin model is proposed for demonstrating cascading failures (break- 
downs, blackouts, collapses, avalanches, ...) that can occur in realistic networks for distribution 
and delivery by suppliers to consumers. A ferromagnetic Hamiltonian with quenched random 
fields results from policies that maximize the gap between demand and delivery. Such policies 
can arise in a competitive market where firms artificially create new demand, or in a solidary 
environment where too high a demand cannot reasonably be met. Network failure in the con- 
text of a policy of solidarity is possible when an initially active state becomes metastable and 
decays to a stable inactive state. We explore the characteristics of the demand and delivery, as 
well as the topological properties, which make the distribution network susceptible of failure. 
An effective temperature is defined, which governs the strength of the activity fiuctuations 
which can induce a collapse. Numerical results, obtained by Monte Carlo simulations of the 
model on (mainly) scale-free networks, are supplemented with analytic mean-field approxi- 
mations to the geometrical random field fiuctuations and the thermal spin fluctuations. The 
role of hubs versus poorly connected nodes in initiating the breakdown of network activity is 
illustrated and related to model parameters. 
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1 Motivation 

Statistical physics studies of cooperative phenomena on random graphs have attracted a lot of new 
attention and undergone impressive new development since it has become clear that many real- 
life interconnected structures are well approximated by random scale-free networks [1, 2, 3, 4, 5]. 
One can say that a paradigm shift is occurring from studies of models for critical phenomena on 
(Bravais) lattices to studies in which such models are defined on random networks. To some extent 
this paradigm shift resembles that from Euclidean geometry to fractal geometry, in the modeling 
of various natural phenomena, but scale-free networks are specific in that they are characterized 
by the presence of a small but important set of hubs. The hubs are highly connected nodes which 
typically have a large influence on the operation and coherence of the structure. 

When we are concerned with the distribution of electricity, or the production and sale of 
material goods by commercial centers to consumers, or the delivery of the daily mail we can 

'Dedicated to Professor Dr. David Sherrington on the occasion of his 70th birthday. 



1 



envisage distribution centers as nodes and consumers as links on a network. The consumers have 
certain demands and the distribution centers certain deliverables. The distribution centers can be 
active or inactive depending on internal conditions and on external criteria such as the demands 
of consumers that are linked to them and the status of other nearby centers. The occurrence of 
hubs in such networks is rather common (main providers, supermarkets, ...). 

We propose a model that can capture three types of common distribution policies in realistic 
trade environments by implementing an appropriate Ising spin Hamiltonian on a random network, 
mainly of scale-free type, but we also consider networks with a scale. We present the three policies 
in a logical order, from obviously cooperative to rather ambivalent, and our main goal in this 
paper will be to demonstrate how a policy in the latter category can cause a breakdown of the 
distribution network. The first policy is one that is guided by demand. The distribution center 
aims at satisfying the consumer demand as closely as possible and will strive to adjust its activity 
to that of neighbouring centers so that the difference between demand and delivery is minimal on 
the links that lead to those centers. Note that every consumer (link) can be served by two centers 
(the adjacent nodes) in this model. This is not unusual. Many consumers rely on an alternative 
option in case their normal provider is not available. A distribution center i is active or up when 
cTj = 1 and inactive or down when ai = —1. The dem,and- guided policy is characteristic of a market 
in which compensation or complementarity is more important than competition. Many physical, 
chemical and biological systems are equipped with similar negative feedback mechanisms, rendering 
operation stable under small enough perturbations. 

The second distribution policy is concerned with a competitive market guided by product quality 
and quantity. In this policy production or distribution centers strive to maximize their activity 
in order to get their deliverables consumed rather than those of their rivals. The centers aim 
at manipulating consumers by creating new demand, preferably in situations where the existing 
demand is low. In this product-guided policy the centers tend to become active especially when their 
neighbours are already satisfying the consumer demand. This policy can therefore be mimicked by 
maximizing the difference between supply and demand, which is opposite to what we postulated 
in the demand-guided policy described before. 

Thirdly, we wish to capture a policy of solidarity rather than competition, in which providers 
display a voluntary or forced shutdown in situations where too high a demand cannot reasonably be 
met. As long as all centers are active there are no problems, but as soon as some become inactive, 
the burden of satisfying the high consumer demand can become too heavy to carry. The centers 
can be unwilling or plainly incapable of rectifying the drop in supply, and increasingly so as time 
progresses and more neighbouring centers fail. This is a pronounced positive feedback mechanism 
which can lead to blackouts in power stations or certain kinds of strikes (when workers are unable 
or unwilling to do more than their peers). A similar conduct may be observed, for example, among 
partners in projects. Here we assume that a link is a project shared between adjacent nodes, which 
can model persons or companies. Companies tend to close or, equivalently, persons tend to stop 
doing their job, when their neighbours no longer invest in a joint project (i.e., when neighbouring 
sites are down). Like in the (second) competitive policy we discussed, in this (third) solidarity 
policy the difference between supply and demand is again maximized, but this time by the supply 
dropping to zero. 

In order to integrate these three policies in an Ising model Hamiltonian description we define, 
besides the spin states associated with the distribution centers, the supply and demand variables 
on the links. For the purpose of this Motivation, we focus on the simplest version of the model, 
but still sufficiently equipped to bring out the essential physics. In later Sections a more refined 
approach will be outlined. For the time being we assume that each distribution center, when 
active, delivers a fixed total amount (set to unity) to all of the links attached, divided equally over 
all links. For the link between nodes i and j with respective degrees Qi and qj, the delivery thus 
equals 




(1) 
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Treating, provisionally, the demand as a uniform constant V throughout the network, we arrive 
according to the policies described above at an energy per link which can be written as 

Eij{ai,aj) = -J[Cij{ui,aj) - Vf, (2) 

where the nearest-neighbour coupling J is negative for a demand- guided policy and posiMve for the 
policies driven by competition or solidarity. The Hamiltonian of the Ising distribution network is 
then the sum of these energies over all links. We obtain the spin Hamiltonian of a given realization 
of a network with A'' nodes, 

<ij> <ij> i=l 

where < ij > denotes nearest-neighbour pairs and the quenched random field Hi acting on the spin 
at node i is given by 

T / 1 \ 

JV (4) 

The spin- independent constant (last term in Eq. (3)) is irrelevant for our purposes and will hence- 
forth be omitted. 

The interpretation of the quenched random field is fairly transparent. For J < (demand- 
guided policy) a distribution center experiences a bias to become active when the total demand of 
the attached links, qiV, exceeds the total average delivery to these links, (1 + J2- q~^)/2, where 
the average is taken over the spin vahies, i.e., over active and inactive states. For J > the spin 
bias is opposite. A distribution center tends to become inactive when the demand is high. 

The nearest-neighbour interaction, J/{2qiqj), is also a quenched random variable, the sign of 
which is equal to that of J. Note that this degree-dependent interaction is reminiscent of that 
in the Special Attention Network model [6], in the sense that high degree is compensated by weak 
interaction. In the distribution model this means that the mutual influence between neighbouring 
centers is proportional to the product of their deliveries to the link that connects them, or to some 
power of this product. We conclude that we are dealing with a random-field (anti-)ferromagnetic 
Ising model on a network for (J < 0) J > 0. The ground state, which minimizes the total energy, 
is in general non-trivial. In addition, "thermal" noise is present, due to occasional maintenance 
shutdowns of centers or fluctuations of their activity caused by other internal factors. We therefore 
consider the model at a flnite temperature T to allow for these more or less random perturbations. 
The ratio ksT / J is a measure of their importance. Preliminary results of Monte Carlo simulations 
of the model at hand were reported by Giuraniuc [7]. 



2 General formulation of the distribution model 

2.1 Construction of a general Hamiltonian 

Consider a static scale-free network with N nodes. The normalized distribution of the degrees, the 
number of links attached to each node, P(q'), with P{q) = 1, is assumed to follow a decreasing 
power-law P{q) a q~^ for large q. The topological exponent 7 usually lies in the regime 2 < 7 < 3 
for real-life networks [3]. Each node i represents a distribution center which can be either active 
(ffi = or inactive (ui = — 1). A link between nodes i and j, denoted by < >, is a consumer 
or a region across which the products are distributed. Contrary to the simplified model proposed 
in the Motivation, the nodes can produce unequal amounts of goods. We assume that the total 
delivery depends on the degree of the supplier. Therefore, the delivery Cij{ai,aj) to a link < ij > 
is defined as 

A,K..) = i.(^i±i) + ^(^), (5) 
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where qi and qj are the degrees of the nodes Unked by < ij >. The exponent jj, controls the total 
production of a supplier. When node i is active or ai ~ +1, supplier i furnishes l/f/f products 
to each of its qi consumers. The total delivery of the active node is thus ql~'^- The case /z = 1 
corresponds to the special case in the Motivation in which all active nodes provide the same total 
supply. Consumers attached to highly-connected nodes then receive a smaller amount. In general, 
the hubs deliver less goods to a single consumer for /U > 0. For = 0, every consumer would 
be receiving the same amount of products in an active state, which also implies that hubs have 
to deliver more goods in total. Finally, ^ < corresponds to the situation in which the highly- 
connected suppliers can deliver more goods to each consumer. 

Following the reasoning introduced in the Motivation, every consumer has a demand V. We 
assume that the demand of link < ij > may depend on q^ and qj . The energy per link is now given 

by 

Eij{oi,aj) = -J[Cij{ai,aj) -Vij] . (6) 
The total Hamiltonian becomes 

N 

^(W) = ^ £;i,(a,,<7,) = -^(i?i+7,(M))<7i (7) 

<ij> i=l 



with 




up to constant (spin-independent) terms, which are neglected. The sums over j run over the 
qi neighbours of node i. This short reasoning leads to two fields applied to each node i. The 
quenched random field Hi is inherent in the network and the interaction field li originates from 
the interactions with spins linked to CTj. 

In this Paper we focus mainly on the third policy introduced in the Motivation. We describe 
avalanches in distribution networks which occur due to the solidarity among suppliers. The in- 
teraction constant J is thus positive. Note that avalanches in anti-ferromagnetic spin systems 
with a uniform external field on complex networks have been studied in Ref. [8]. Avalanches in 
distribution models were also studied using sand piles on networks [9] and in various models based 
on critical loads [10, 11, 12, 13]. Recently, also models studying percolation in interdependent 
networks were used [14]. 




2.2 The demand function 

Finally, we propose a suitable demand function. The demand of a link < ij > between nodes i 
and j is a combination of a link-dependent and a homogeneous part: 

(8) 

where a(> 0) and b are real constants. The notation (•) denotes an average over the degree 
distribution P{q)- The first term, the supply- adjusted demand is a consequence of the form of the 
fields Hi in Eq. (7). The parameter a controls which fraction of the normal capacity is demanded by 
the consumers on the link. Each link < ij > also features an intrinsic uniform demand, independent 
of the deliveries of the nodes i and j. The second term, regulated by the constant 6, provides such 
a global demand. The onset of a homogeneous consumer demand [b ^ 0) results in a total field 
on the network which has the same order of magnitude as the interaction field in an active state 
averaged over an ensemble of different networks. Note that in certain systems individual consumers 
can deliver goods to the market. For instance, individual households can contribute to the power 
grid with the output of their own solar cells. Negative values of b are therefore not excluded. 
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The total field acting on node i is now given by 



3 = 1 



(9) 



Note that only the last term depends on the interaction with the spins on the neighbouring sites. 
It will play a prominent role in the cascading effect. For general values of /i, both the interaction 
and the quenched random field applied on node i will depend on the connectivity of node i and on 
the connectivity of its neighbours. However, in some cases simplifications can be made. If = 0, 
the parameters a and b can be replaced by a single parameter ao = a + 6/2. The total field is then 
simplified to 

J 

Hi+h{{a}) = Jqi{l-ao) + -y2(Tj; for /x = 0. (10) 

^ i=i 

Both the interaction field and the quenched random field acting on a node are then independent 
of the degrees of the neighbouring sites. For arbitrary fi the quenched random field is indepen- 
dent of the connectivity of the neighbouring nodes for a = 1. Note finally that for a = 1 and 
6 = the demand just matches the average supply (the average being taken over active and in- 
active providers) and the quenched random field vanishes, i.e.. Hi = 0, Vi. An Ising model with 
connectivity-dependent coupling constants in the absence of an external field was studied by Giu- 
raniuc et al. [15]. It was shown to be equivalent to a model with homogeneous couplings and with 
a modified topological exponent 7' = (7 — /i)/(l — /i). 

The quenched random and interaction fields featured in Eq. (9) depend on the local structure of 
each network. Handy analytic approximations can be obtained by performing an ensemble average 
over different network realizations with the same degree distribution. In a first step, we average 
over the degree of the neighbours to obtain the average field applied on node i. Let Pn^ (qj) denote 
the probability that the neighbour of node i has connectivity qj . We assume that the network is 

uncorrelated. The distribution Pn\qj) is then independent of node i and related to the degree 
distribution P{qj) by the relation Pn{qj) = <ljP{Qj) /{<!)■ Therefore we define the average over the 
degrees of the neighbouring sites of node i as follows: 

«-»^ = flE^- (11) 

This type of averaging is a "topological" mean-field approximation. Simultaneously, we may in 
some situations wish to apply a "thermal" mean-field approach to the thermally fluctuating spin 
variables and replace the actual spins of the neighbours by a mean spin. If we furthermore assume 
that the mean spin is homogeneous throughout the network, the average quenched random and 
interaction fields acting on node i in the mean-field approach are 

= Jq]-^^l^ + Jql-^^^(l-a-h), (12a) 

«/i({(7})»i = Jql-i'^^^aa,, (12b) 

with Gav tfie average spin of a node. Deviations from these averages are important for a network 
in which the degree distribution P{q) possesses diverging moments. For instance, in a scale-free 
network with topological exponent 7 < 3, the second moment diverges. Note that the average 
quenched random and interaction fields fields on a node diverge for < 2 — 7. Therefore we will 
henceforth limit our attention to the case > 2 — 7. Since we only consider scale-free networks 
with a finite mean degree, i.e., with 7 > 2, this limitation is only important for /i < 0. 
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In a second step, we average over the degrees of the nodes using the degree distribution P{q). 
The average Hamiltonian (^{uav)) is then given by 

{n{aav)) = Y.Y. Pi(liy<^v ^Hi + Ii{aav) >i • (13) 

i Qi 

Inserting both members of Eq. (12) we obtain 

{H{aav)) = -NJaav ((a^"'^)^ + -^^^ (1 - a - 6 + a^,)) . (14) 

Eq. (14) provides a mean-field expression for the Hamiltonian averaged over different reahzations 
of the same network. 



3 Requirements and favourable circumstances for a network 
collapse 

For the study of a collapsing network, driven by competition or solidarity (J > 0), we exploit 
the ferromagnetic phase which appears in the Ising model at sufficiently low temperatures and for 
sufficiently weak random-field fluctuations. Recall that the ground state (T = 0) in the absence of 
fields consists of two degenerate configurations: the active network in which W i, ai = +1, and the 
inactive network in which V i, CTj = —1. By applying a uniform bulk field, the degeneracy is lifted 
and a unique equilibrium state emerges. Nevertheless, the network can reside for a long time in 
the oppositely magnetized "mctastable" state if the fields and (thermal) spin fluctuations are too 
small to trigger the collapse to the equilibrium state. 

Metastability is defined dynamically in our context as partial stability with respect to certain 
perturbations, in contrast with absolute stability (with respect to any perturbation). In our model, 
in the presence of random fields of microscopic origin, which are inhomogeneous and regulated by 
the parameters a, b and fj,, the ground state may be non-trivial and will coincide with the all- 
up or all-down states only under certain conditions on the random fields. Therefore, considering 
metastable states, we need to distinguish between states that decay to a ferromagnetic state or 
evolve to a "glassy" state. We will be interested mostly in the ferromagnetic state because in a 
glassy state the network is still more or less active. 

Thus, in our study of distribution networks driven by solidarity, we will mainly focus on the 
decay of the all-up state (active network) to an essentially all-down state (inactive network). When 
the global consumer demand (modeled by b) is increased, the metastable states we study by means 
of a suitable spin-flip dynamics, mimic the behavior of certain realistic systems. For example, from 
1988 to 1998, US electricity demands increased by nearly 30% while the network capacity grew 
only by 15% [16]. Apparently as a result of this widening gap between demand and supply, the 
system became metastable, which only became apparent when a large part of the power grid broke 
down. 

In the following, we start from an active network and determine under which conditions a 
collapse to the inactive ground state is likely to occur. After an active period, avalanches to the 
inactive state can only take place if the system initially resides in a metastable active state and 
then decays to the energetically more favourable inactive state. In the following these requirements 
are converted into conditions in terms of the constants of the distribution model. 

There are indeed two basic restrictions in order to sec avalanche effects. The first one concerns 
the metastability of the active state. The active state should remain intact for a siifficiently long 
time. Focussing on singlc-spin-flip dynamics, this requirement corresponds to the impossibility 
that at zero tempcratm-e a spin spontaneously flips from -1-1 to —1 while its surrounding remains 
active (-1-1). In such a hypothetical process, only the local energy associated with node i changes 
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by an amount AE?-^ given by 

AE'/ =2{Hi + Ii{aav==+l)). (15) 

At zero temperature, single-spin flips are excluded provided AE^^ > 0. The definition of the fields, 
Eq. (9), then leads to the requirement 

(q) I 1—a 2 — a 1 I ,,. „N 
b < , ^ —n- + > -u , Vz (16) 

where the sum is over the Qi neighbours of node i. The condition depends on the specific local 
structure around each node in the network. 

In order to obtain a simple and useful analytical approximation to this, we can average over 
the degree of the nearest neighbours. The averaging procedure defined in Eq. (11) entails the 
replacement 

1 ^ (gi-'^) 



r < (9) 



(17) 



Note that the substitution defined in Eq. (17) is only exact for /x = 0. For /z ^ it is a topological 
mean- field approximation. Within this approximation the active state is metastable (i.e., stable 
against single-spin flips at T = 0) if for all possible degrees m, 



(g) (l-g) 



b< .T-,,/ , ' +2-a,yi. (18) 



This condition is equivalent to 

<^Hi + Ii{aav = +1) >i > 0, Vi, (19) 

which is the metastability criterion for the mean fields defined in Eq. (12). The mean-field approx- 
imation can therefore by formulated directly in terms of the mean fields. 

Eq. (18) clearly sets an upper limit to the global demand. If it is too large, a metastable 
active state is not possible and there will be an "immediate" decay to the inactive ground state. 
Depending on the signs of 1 — a and /i, two situations can be distinguished. For /i > and a < 1, 
it suffices that the largest possible degree, K, satisfies Eq. (18). On the other hand, for /u > and 
a > 1, the smallest possible degree, m, is the important one. For /x < 0, the converse is true. 

If more than one spin is allowed to flip at the same time, the criterion Eq. (18) must be extended. 
The energy difference associated with a multiple-spin-flip process can be smaller in magnitude than 
the energy difference pertaining to a single-spin-flip process. For instance, for a process in which 
two nearest neighbours i and j flip, the total energy difference of the double-spin-flip process AEfj 
would be 

AEf^ = AEf + AEf - 2J(g,g,)-^ (20) 

which (for J > 0) is not only always lower than the sum of the energy differences of the two 
single-spin-flip processes separately {AE^-^ + AE^^), but can even be lower than the energy change 
involved in a single-spin-flip process. This can be appreciated by considering for example the 
special case a = 6 = 1, for which the average value 4C AE^^ vanishes. The energy difference 
may in general be reduced further if more than two neighbouring spins can flip simultaneously. In 
real distribution networks multiple-spin flips can model multiple suppliers failing at the same time, 
when, for instance, raw materials are exhausted. However, we will not consider these instances but 
rather allow for random technical failures besides deliberate decisions arising from interactions, 
making it unlikely that multiple suppliers break down exactly simultaneously. Therefore, for 
simplicity, single-spin-flip dynamics is assumed throughout the remainder of our paper. 
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As a second restriction, we require the inactive state to have a lower total energy at T = 
than any other state. This condition ensures that the breakdown occurs towards the inactive 
state rather than to another, glassy state, which minimizes the energy in regions of the parameter 
space where the local random fields dominate the ferromagnetic interaction. The reason for this 
limitation is that we wish to study principally the blackout phenomenon in which, after some time, 
a state with very low activity is reached. 

A sufficient condition for the inactive state (all spins down) to be the ground state (at T = 0) 
is that all local quenched random fields have the same sign, in this case negative. The exact 
requirement is ffj < 0, Vz, or 



^>(l-«)7I^hF + TE7Fh (21) 




In the mean-field approximation adopted in the previous paragraph this becomes 



-Ol^^ + l-". «. (22) 



Eq. (21) and its simplified version Eq. (22) provide a threshold for the demand sufficient to observe 
decays to the inactive state. When the demand is smaller than this threshold, the absolute stability 
of the inactive state is not guaranteed. Then the system may remain stable in the active config- 
uration (all spins up) or evolve to a glassy state. Once more, different regimes can be identified 
according to the vahies of a and /n. For example, for positive /i, if a < 1 it suffices that the smallest 
possible degree satisfies Eq. (22). However, if a > 1, the relevant quantity is the largest degree. 
These statements are to be interchanged if ji is negative. 

The two conditions Eq. (16) and Eq. (21), or their mean-field approximations Eq. (18) and 
Eq. (22), determine the region in a (a, 6)-phase diagram where collapses of an active network to an 
inactive one should be observable. A graphical representation of the phase diagram for distribution 
models with fj, — and ii — 0.2 can be found in Fig. 1. The shaded region on the Figure indicates 
the ranges of a and b within which the two requirements are satisfied for a scale-free network with 
topological constant 7 = 3, minimal degree m = 2 and 1000 nodes. The uppermost line (dotted) is 
the numerically exact upper bound for the stability of the active state against single-spin flips at 
T = 0, the so-called metastability limit, given by Eq. (16). This exact condition was determined 
by simulation of the model on scale- free networks. 

The networks were generated in two steps using the uncorrelated configuration model. In a 
first step, A'^ nodes are created, each with their own degree, chosen randomly according to the 
distribution P{q), which is a decreasing power law. The minimal degree is m = 2, while the 
maximal degree is set to K = ^/N. For 7 > 3 this is not a restriction, since the maximal degree 
satisfies K oc A^i/(t-i) [17]. On the other hand, for 7 < 3 this restriction avoids correlations in the 
network [18]. In a second step, links are laid randomly between all the nodes, with the constraint 
that at the end of the linking procedure every node should have the degree it was given in the first 
step. Self linking and multiple linking are avoided. The results obtained from different network 
realizations differ slightly. In practice, averaging over 10 networks is sufficient to obtain accurate 
reproducible results. In this manner we determine the values of a and b for which Eq. (16) and 
Eq. (21) are satisfied for all nodes. The latter condition leads to the second uppermost line (solid) 
in the figures (straight line for /x = 0; line consisting of two straight segments, with a break at 
a = 1, for /J, ^ 0). Above this line the inactive state is, with certainty, the groimd state (at T = 0). 
The region in between the two lines discussed so far is susceptible of avalanches of spin flips, i.e., 
blackouts of the network, and is shown shaded (in dark grey). 

The model possesses an overall symmetry, which is reflected in the phase diagrams in Fig.l. 
Inspecting Eq. (7) and Eq. (9), we conclude that the full Hamiltonian is invariant under the 
transformation: a — )• a' = 2 — a; 6 — )• 6' = —b; cr — )• cr' = —a (flipping all the spins). This 
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symmetry implies that the lines in the phase diagram pertaining to the active state (all spins up) 
can easily be drawn by applying the above transformation of b and a on the lines associated with 
properties of the inactive state (all spins down). In particular, the lowermost line (dashed) marks 
the metastability limit of the inactive state, and the second lowermost line (dash-dotted; broken 
for fi 0) indicates the limit of absolute stability of the active state. Below this line the active 
state is, with certainty, the ground state (at T = 0). A remarkable feature of the phase diagram 
now emerges. While for = the two limits of absolute stability coincide, so that there are only 
(two) ferromagnetic ground states, a new zone appears for /i 7^ 0. In this zone, which wc term 
"no-man's land" , the ground state is not with certainty ferromagnetic. It may be a glassy state, 
characterized by local random fields of either sign (up or down) that try to orient the spins along 
their direction, at low T, as they compete with the ferromagnetic spin-spin coupling. The no-man's 
land is shaded in light grey. Note that the width of this no-man's land vanishes for b = 0, and 
also, for a = 1, for which the random fields are trivially all of the same sign, determined by the 
value of the remaining free parameter a or 6, respectively. 

In general, the phase diagram depends on the network topology. However, the case /x = 
provides an exception. As mentioned before, the single parameter ao = a+b/2 suffices for describing 
the supply and demand functions, if /i vanishes. The conditions Eq. (18) and Eq. (22), which are 
exact for /x = 0, then correspond to the simple conditions 1 < ao < 3/2. Note that these can also 
be obtained directly starting from Eq. (10). These conditions imply that there are no distinctions 
between different networks at zero temperature if /i vanishes. Indeed, for /i = 0, the sign of 
the quenched random field is independent of the network structure, as can be seen in Eq. (10). 
Therefore, the topology of the network does not affect the state of the system at T = 0. This 
prediction is verified by comparing simulations on scale-free networks with simulations on random 
Erdos-Renyi graphs. In an Erdos-Renyi network all pairs of edges have equal probability to be 
linked, which implies a Poissonian degree distribution. We verified that the condition Eq. (16) 
leads to one and the same straight line in the (a, 6)-plane on an Erdos-Renyi network and on a 
scale-free network with 7 = 3. The same is true for Eq. (21). However, if 7^ 0, the network 
topology exerts a crucial influence on the phase diagram. We verified that, for example, both with 
respect to the mean-field conditions and the exact ones, the phase diagram is different for networks 
with different values of 7 [19]. 

In the following we go into more detail as regards the interesting comparison of the mean-fleld 
approximation and the (numerically) exact results for the boundaries in the phase diagram. As we 
already noted, the topological mean-field approximation is (only) exact for n = 0. The simulations 
confirm this. To estimate the quantitative difference between, for example Eq. (16) and Eq. (18) 
we calculate the variance, over all network realizations, of a quantity associated with the sum in 
the last term of Eq. (16). We obtain 



This result allows us to estimate the importance of the random field fluctuations. Clearly, the hubs 
are least affected by the topological disorder since the SD scales with their degree qi as l/^/qi. The 
amplitude of this power law depends on /z, in a manner which is easy to interpret. For example, 
for fj, = —1,0.5 or 1 the variance is proportional to {q){q'^) — (9^)^, (q) — (9^^^)^ or (q^^) — {q)~^, 
respectively. 

We now proceed to assess quantitatively the difference between the exact metastability limit 
Eq. (16) and the approximate one Eq. (18) for /x = 0.2. Clearly, as Fig. 2(a) shows, the mean-field 




(23) 



Consequently, the ratio of the standard deviation to the mean is given by 




(24) 
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approximation leads to a continuous piece- wise linear curve which is broken at o = 1, since for 

a < 1(> 1) the right-hand-side of Eq. (18) is minimized by the maximal (minimal) degree present 
in the network. Interestingly, the exact (dotted) curve for /x = 0.2, shown in Fig. 2(a) and also in 
Fig. 1(b), is a straight line, without any singularity, since the second term in the right-hand-side 
of Eq. (16) is typically minimized by a node i of degree (jj = 2 connecting two hubs {cjj ^ 1). 
This minimal value lies some 3 to 4 standard deviations below the mean. Consequently, the entire 
right-hand-side of Eq. (16) is minimized by one and the same node i with a low value of qi for 
all a G [0,2]. For /i larger than some threshold this is no longer the ease and the numerically 
exact metastability limit is no longer a straight line but a gently bent concave curve. Still, its 
shape appears smooth and is therefore qualitatively different from the broken curve found in the 
mean-field approximation. This is conspicuous in Fig. 2(b) where both are shown for /i = 1. 

We repeat our analysis for the absolute stability limit of the active state. This is the curve below 
which the active state is with certainty the ground state (at T = 0). It is derived numerically exactly 
by applying the transformation a — >■ a' = 2 — a; 6 — )■ 6' = — 6 to Eq. (21) and in the mean- field 
approximation by applying the same symmetry to Eq. (22) , which in both cases simply reverses the 
inequalities concerned. In contrast with the metastability limit, the slope of the absolute stability 
limit displays a discontinuity at a = 1 for both the exact and the mean- field versions, as can be 
seen from the fact that both Eq. (21) and Eq. (22) contain the prefactor 1 — a. For a < 1 the 
limit is defined through a hub (typically the node with the highest degree) and for a > 1 through a 
node with a low degree. For the former case, the mean-field approximation is accurate as expected 
since for hubs random field fluctuations are small, being proportional to 1/y^. For the latter case, 
random field fluctuations are more important and, indeed, a clear difference emerges between the 
mean-field upper bound and the exact one. In Fig. 2(a) we compare the exact and the mean-field 
curves for /j. = 0.2, both for the metastability limit of the active state and for the absolute stability 
limit of the same state. 

The region in which avalanches can occur depends sensitively on the exponent fi. Comparing 
Fig. 1(a), Fig. 1(b) and Fig. 2(b), the size of the region in which blackouts are possible shrinks 
as increases. The same trend is observed when /x is decreased from zero [19]. We conclude 
that the region in parameter space in which our two criteria for avalanches are fulfilled shrinks 
as |/Lt| increases. This region is also sensitive to the topological exponent 7. For instance, for 
/X = 1, when 7 is decreased from 5 to 2 there appear typically more nodes with higher degrees. 
Consequently, the mean degree (q) increases. This raises the absolute stability limit of the inactive 
network for a < 1 and lowers the metastability limit of the active state for a > 1, while leaving 
the other segments of the phase boundaries unchanged (as can be seen qualitatively by inspecting 
the equations for the mean- field (meta-)stability criteria). As a consequence the region susceptible 
of avalanches is squeezed more tight about the center (a = 1,& = 0) of the phase diagram. At 
this point we would like to mention other studies in which the resistance against cascades was 
studied as a function of geometrical disorder. In models based on critical loads and percolation in 
interdependent networks it was found that a heterogeneous network is less resistant against higher 
loads than a homogeneous one [12, 13, 14]. In our model degree heterogeneity (low 7) appears to 
strengthen the network at low a but to weaken it at high a. 

Also the size of the network influences the region available for collapses, mainly through the 
maximal degree which depends on N. The region in which avalanches are observed decreases if 
more suppliers are present in the network, as can be seen from Eqs. (18) and (22). However, a 
numerical evaluation of these criteria indicates that the effect of the number of suppliers is rather 
small compared to the effects of the values of fi and 7 [19]. We conclude that the ranges of the 
demand parameters a and b for which collapses occur are influenced most strongly by the network 
topology and the degree dependence of the delivery in the distribution system. 
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4 Collapse properties 



4.1 Distribution model at finite temperature 

In real-life distribution systems, individual suppliers can fail to deliver their goods, for instance due 
to a defect or malfunction in the production process. After repair the delivery resumes. Therefore, 
wc introduce a generic "temperature" T in the distribution model to quantify the rate of random 
spin fluctuations, from the active to the inactive state and back. At sufRciently high temperatures, 
the mean spin magnetization M = J^i'^i/^' which represents the mean network activity, tends 
to zero as a function of T according to a Curie- Wciss-typc law, i.e., M cx 1/T for T oc' (cf. a 
paramagnet in a small external field). We will, however, focus mainly on lower temperatures, still 
in the ferromagnetic regime, for which a network failure can occur. 

Wc perform simulations on scalc-frcc graphs which arc constructed using the iincorrclatcd 
configuration model introduced in the previous Section. All networks contain 1000 nodes, the 
minimal degree of which is m = 2 and the maximal one K = \/T000 « 32. Spins are updated 
using singlc-spin-flip dynamics with the Metropolis updating rule. The simulations start from a 
metastable active state. Collapses can be found for various values of the parameters n, a, b of 
the distribution model and for different values of the topological exponent 7. Some examples are 
shown in Fig. 3, in which the time evolution of the magnetization is plotted. A single time step 
corresponds to 1000 single-spin fiips. Different types of collapses are observed in different ranges 
of model parameters, as we will illustrate in the following. 

A first point of attention concerns the nodes which initiate the collapse. Apart from the mean 
magnetization of the network. Fig. 3 also shows the mean magnetization of the highly connected 
nodes (with at least 8 links) and that of the poorly connected nodes (with only two neighbours). 
In the first and second collapses. Fig. 3(a) and Fig. 3(b), the poorly connected nodes exhibit the 
largest fluctuations before the collapse. The nodes with a low degree thus initiate the breakdown. 
However, the actual collapse only takes place when also the hubs start to flip. As long as the hubs 
remain active, their large influence in the network prevents a blackout. The opposite behavior is 
found in the avalanches shown in Fig. 3(c) and Fig. 3(d), in which the hubs initiate the collapse. 
Again, the mean magnetization remains positive until also the multitude of less-connected nodes 
starts to collapse. Thus in all four cases only a certain subset of nodes initiates the network 
collapse, but the network undergoes a full breakdown as a consequence of a collective effect, in 
which all nodes are involved. Our model thus displays the cooperative character of distribution 
networks as described in the Motivation. 

A simple mean-field argument suggests a condition for determining whether the hubs or the 
poorly connected nodes initiate the collapse. The collapse is initiated by the nodes which have the 
largest fluctuations in the active state. We introduce an activation temperature for a node with 
degree q, Tactiq), by equating its thermal fluctuation energy to the average energy needed to spin 
flip a node with degree q from -|-1 to —1. When T > Tactio) the spin of the node will undergo 
significant thermal fluctuations. Within the mean-fleld approximation, Tact{(l) is given by 

kBTacMi) = 2 < iJ, + h{{a}) »i= Jq]-^''{1 - a) + Jq]-" (1 - a - 6 + aav) ■ (25) 

Wc focus on the degree-dependent activation temperature when the network is in an (almost) 
fully active state, thus with aav ^ 1- Then, the nodes with the smallest Tact (9) display the 
largest fluctuations in the active state and can thus initiate the collapse of the network when T 
is slowly raised from zero. Note that whether the activation temperature is either increasing or 
decreasing as a function of the node degree, depends on the signs ofl — o, /U— 1/2 and 2 — a — 6 
(assuming aav = 1)- Let's test these ideas against some of the simulations of network breakdowns 
shown in Fig. 3. For example, using the model parameters associated with Fig. 3(a), wc obtain, 
with aav ^ 1, ksTactiQ) ~ J(0.30g'''^ -|- 0.31 g°-^), which implies that the activation temperature 
increases with the node degree. Therefore, the poorly connected nodes should initiate the collapse. 
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which is confirmed by the simulations. Taking the parameters associated with Fig. 3(d), we 

find ksTactiq) ~ J(0.28 + 0.15 g"^), so that the highly connected nodes should show the largest 
fluctuations in the active state. The simulations confirm that the hubs indeed initiate the network 
collapse. 

Information about the type of nodes that initiates a collapse is of great interest in real-life 
networks. If hubs tend to be the most fragile and thus most fluctuating suppliers, it is best to 
invest more effort in protecting them rather than the least connected nodes. Of course the converse 
is true if the poorly-connected nodes are more prone to initiate the collapse. The criterion of 
Eq. (25) could therefore be used to strengthen networks purposefully against the consequences of 
accidental malfunction of distribution centers. 



4.2 Effective strength of thermal fluctuations 

As can be seen in Fig. 3, the magnitude of the thermal fluctuations depends not only on the 
value of T but also on that of other parameters of the distribution model. Even if the demand 
is fixed (constant a and 6), thermal fluctuations may still be reduced or amplified depending on 
the network structure, through the topological exponent 7, and depending on the delivery in the 
model, through the delivery exponent 11. The starting point of our further analysis is the mean 
Hamiltonian, Eq. (14). The expression for this Hamiltonian together with the expressions for the 
amplitude of the fluctuations of the non-local term in the quenched random fields, Eqs. (23) and 
(24), suggest that, subject to conditions to be specified, the dependence on the network topology 
and on fx might be captured by the single parameter {q^~^)'^ / {q)- This prompts us to test the 
conjecture that the mean Hamiltonian may possess the following scaling property, 

{n{aav))^-NJaav^^l^f{.aA<yav), (26) 

where the function /(a, 6, (Jav) is independent of jj, and 7. This Ansatz is most likely to be valid 

in at least one of the two following circumstances. Firstly, if a is sufficiently close to one (a ss 
1), the first term in Eq. (14) is negligible compared to the second one and Eq. (26) holds with 
f{a,b,<Tav) ~ — 6 + <Tav The second situation occurs for large networks and small enough /x, i.e., 
< 7 - 2, which leads to (g^"^^) ~ {q'^-'')'^ / (q) , as can readily be shown analytically. Indeed, 
in the thermodynamic limit, N,K — >• 00, and converting sums over the degree distribution to 
integrals [20], one finds 

Note that, in view of Eq. (24), the condition <C 7 — 2 thus ensures that the random-field 
fluctuations are small. Under these circumstances Eq. (14) indeed takes the scaling form of Eq. (26) 
with /(a, 6, ffav) = 2 — 2a — 6 + Uav Although, strictly speaking, this last result is only valid for 
the range of /z specified above, numerical inspection shows that it is a good approximation even for 
values of /i up till about unity, provided the network is large enough for the /T-dependence of the 
averages to be negligible. In general, the approximation remains useful also for finite networks, as 
numerical analysis shows. In the remainder of this section we assume /x > 0. 

According to Boltzmann statistics and using the Ansatz of Eq. (26), the probability to observe 
a network with mean spin aav satisfies 



f NJ{q^-'^)^a,J{a,b,aa^) 

^^"'^^^ {-T^m 2 — 



(28) 



Apart from the constants related to the demand function, all model parameters {n, 7 and T) are 
only present in the first factor in the exponential function. We therefore absorb the dependence 
on the delivery exponent, the topology and the temperature into a single parameter 6, defined as 

e = TT^r- (29) 
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For systems with fixed a and 6, the finite-temperature behavior of the system is thus controlled 

by 0, which acts as an effective tem,peratiire. In the next subsections, the; cfFcicts of both the 
delivery exponent fjt and the topological exponent 7 on finite-temperature collapses are investigated 
separately. 

4.2.1 Effect of the delivery exponent /i 

We now focus on a distribution model with fixed demand constants a and 6 on a scale-free network 
with fixed topological constant 7 and investigate the effect of different values of the delivery 
exponent n G [0, 1]. Since {q^~^) decreases with increasing /i, the effective; tciniperature of the 
network defined in Eq. (29) increases as a function of fx. Thermal spin fluctuations are thus larger 
in a distribution model with larger /z than in a model with smaller /z at the same temperature T. As 
/i increases, the network tlms becomes more vulnerable, since the collapse temperature decreases. 
The considered distribution network could therefore be effectively strengthened by a decrease of ji, 
i.e., by rendering the amounts of goods delivered to each customer more homogeneously distributed. 

In addition to this "mean-field" effect there is an effect of the fluctuations of the quenched 
random nearest-neighbour coTiplings, or random-bond fluctuations. These are of topological origin 
and induced by the (quenched) degree fluctuations of the nodes. To understand this we recall that 
previous studies have shown that the (equilibrium) critical temperature Tc of the model in zero 
external fleld depends rather sensitively on the values of 7 and /x. In particular, for /x = 0, the 
critical temperature is finite as long as the second moment of the degree distribution is finite, but 
diverges as a function of the network size N when {q'^) diverges [21, 4], which is the case for 7 < 3. 
For ji ^ Q, whether or not Tc is finite (for an infinite network) is determined by whether or not 
the effective exponent 7' = (7 — Ai)/(1 ~ m) exceeds the value 3 [15]. Therefore, it is important 
to check also the value of 7' in our distribution systems. In particular, when /i is increased (from 
towards 1), the exponent 7' increases and the hubs become less numerous and less pronounced. 
Consequently, Tc decreases and this also renders the thermal spin fiuctuations more important. 
This effect will be most relevant for networks with 7 < 3 and 7' > 3. 

The effects we described are confirmed by simulations in which the mean magnetization is 
studied versus temperature. In each such simulation we initialize our network in the metastable 
active state and we update the spins using single-spin flips during 4000 time steps. At sufficiently 
low temperatures, the metastable state remains stable during the entire simulation and the final 
magnetization (after 4000 time steps) remains close to one. Repeating this procedure for a sequence 
of fixed temperatures, upon increase of the temperature one observes a transition to the regime 
in which the active state is no longer metastable on the time scale of the simulations. This 
(non-equilibrium) transition, which takes place at a breakdown temperature Tj, is marked by a 
magnetization jump to a negative value. An inactive state, with Uav ^ —1, is then reached before 
the end of the simulations. At still higher temperatures, the final magnetization becomes smaller 
in absolute value and displays a Curie- Weiss behavior reminiscent of the paramagnetic state. Note 
that, in the absence of symmetry-breaking fields, the final magnetization would approach zero 
rather sharply at the equilibrium critical temperature T^. We also verified this in our simulations. 
Obviously, T^kT^. 

Simulation results for different values of ji are shown in Fig. 4(a), for a = 1, and in Fig. 4(c) 
for a ^ 1. The magnetization curves tend to be stretched and shifted as pi decreases, reflecting the 
fact that networks remain stable up to a higher temperature for smaller values of ji. Interestingly, 
for a = 1, if we plot the magnetization as a function of the effective temperature 0, a good data 
collapse occurs, as is conspicuous in Fig. 4(b). Not only do the different curves for T > T;, in 
Figs. 4(a) fall onto a single curve in Fig. 4(b), but also the different values of lead to practically 
one and the same value of 8^. We argue that the high quality of the data collapse has two reasons. 
Firstly, for a = 1 the scaling Ansatz Eq. (26) is properly valid (at mean-field level) and secondly, 
the equilibrium critical temperatures for the different networks (in zero field) are not very different. 
(An exception could, in principle, have occurred in the borderline case 7 = 3 and /U = 0, for which 
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To diverges. But this divergence is very slow, in the manner Tc oc logN, for large N.) Notice how 
the quahty of the data collapse degrades if a 7^ 1, as can be seen by comparing Figs. 4(b) and 4(d). 
In addition, there is now also a significant spread on the values of 6(,. These effects could have 
been expected, since the scaling Ansatz Eq. (26) is not well satisfied for a 7^ 1, if is arbitrary. 

Wc conchidc that, at least in the range < /i < 1, the effects of varying and varying T are 
not independent. A change in jj can be compensated or "absorbed" by a change in T. This is most 
clearly so for a w 1 and for 7 sufficiently large (i.e., in practice for 7 > 3). 

4.2.2 Effect of the topological exponent 7 

In a second application of the use of the effective temperature of Eq. (29), we determine the 
effect of the network topology on the network activity or spin magnetization. We focus on the 
distribution model on scale- free networks with m = 2, 1000 nodes, constant a, b and /x and we 
vary the topological exponent 7. Converting sums over the distribution function into integrals [20] , 
Eq. (29) leads to 

Q_ (7-2 + M)^ (K^-^-m^--)(A-i-^-mi-'^) ^ 

(7-l)(7-2) (^-2-7-^ _ j„2-7-/.)2 

As in the previous subsection, we distinguish between effects at mean-field level and effects of 
the quenched random-bond fluctuations on the characteristic temperatures T{, (breakdown) and 
(bulk criticality). 

At mean-field level, the behavior of O as a function of 7 is complex and depends both on the 
parameter /x and on the size of the network through the maximal degree K. However, when = 
or /X > 1, some simplifications apply. For = 0, oc l/(g) so that 6 increases with increasing 7, 
regardless of m and K. For vanishing n, the network thus becomes more vulnerable to collapses 
if 7 increases, i.e., when there are less hubs in the network. In such a regime, hubs produce more 
goods than poorly connected nodes. The network thus appears more resilient against collapses 
when there are more hubs and when they are more productive. 

The opposite behavior occurs for /x > 1. For constant temperature T, Q decreases with increas- 
ing 7 and thus the network becomes more vulnerable to collapses when 7 is decreased. The system 
is thus less prone to failure when there are less hubs and more nodes with a small degree. For large 
/i, nodes with small degrees provide the larger quantities of goods. Between the two regimes, i.e., 
for < /X < 1, there is a complex transition region in which the behavior of as a function of 7 
depends more subtly on the value of /x and the size of the network. Numerical simulations indicate 
that the effect of topology on the thermal fluctuations is rath(;r small in this regime. 

The effect of quenched random-bond disorder has already been discussed in the previous subsec- 
tion. It suffices to recall that the value of the effective topological exponent 7' determines whether 
wc are dealing with a strongly expanded temperature scale {Tc — >■ 00) or a normal one (with finite 
Tc). 

We illustrate the above qualitative features with simulation results in Fig. 5. In Fig. 5(a), the 
simulations use a distribution model with n = 0, while in Fig. 5(c) /x = 1 is taken. In both cases 
a = 1 is assumed so that the mean-field Hamiltonian has the simple scaling form of Eq. (26). For 
/X = network collapses occur at lower temperatures for networks with larger 7, which implies 
that networks with more hubs and more highly connected hubs are more robust against thermal 
fluctuations. The opposite behavior is observed for /x = 1 (see Fig. 5(c)). The simulations thus 
confirm our expectations based on the dependence of © on /x and 7. In Figs. 5(b) and 5(d), the 
magnetization of the networks for the systems of Figs. 5(a) and 5(c) is plotted as a function of the 
parameter 0. 

In Fig. 5(b) the data collapse is far from perfect. This is at first sight surprising because the 
conditions a = 1 and /x = seem ideal prerequisites from the point of view of the validity of the 
Ansatz of Eq. (26). However, the networks examined are qualitatively different in the sense that 
for 7 < 3 the degree fiuctuations are important ((g^) diverges for N — >• 00), while for 7 > 3 they 
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are not. If we take into account the finite network size, we obtain the following estimates for the 

equilibrium critical temperatures in zero external field, using the analytic results of previous works 
[15, 4]. For 7 = 5, 3, and 2.2 we find fc^Tc/J « 0.60, 1.95 and 3.89, respectively, which spans a 
broad range. These values appear consistent with the behaviour of M(T) (in non-zero external 
field) shown in Fig. 5(a). The simple "mean-field" scaling underlying the definition of the effective 
temperature Q is not quite sufficient to suppress the rather large effect of the degree fiuctuations 
for low 7. This is why the magnetization curves and the breakdown temperatures show only a 
rather poor data collapse in Fig. 5(b). 

In contrast, a much better "universality" is clearly emerging in Fig. 5(d), which is for systems 
with /It = 1. For these systems 7' = 00 so that the networks all behave effectively as Poissonian 
networks, with a finite Tc which scales in a simple manner with (q). It is conspicuous in Fig. 5(d) 
that both the M{T) and the effective breakdown temperatures coincide well for all three cases. 

5 Conclusions 

In this Paper, we introduced a model to describe cooperative behavior in various real-life distribu- 
tion systems. We specialized to networks in which a policy of competition (attempting to create 
new demand) or solidarity (reluctance to meet too high a demand) among suppliers takes effect so 
that there is a tendency towards maximizing the gap between delivery and demand. The resulting 
positive feedback mechanism can cause an initially active network to function during a certain time 
and then collapse to an inactive state. To implement this our model utilizes an Ising-spin system 
with quenched random fields. The couplings are ferromagnetic in order to describe the positive 
feedback policies: spins align preferentially with their neighbours (up: active state; down: inactive 
state). Each node, of degree q, models a supplier which has a certain degree-dependent delivery 
controlled by an exponent fi and each link models (a number of) consumers with a certain demand, 
controlled by the amplitudes a (supplier-adjusted demand) en h (global demand). The degree dis- 
tribution P{q) is (mostly) of power-law type, pertaining to a scale-free network with topological 
exponent 7. The system is studied analytically, in part by using mean-field approximations, and 
numerically using Monte-Carlo simulations with the Metropolis updating rule. 

Avalanches between an active and an inactive state are typically only observed if the inactive 
state has a lower total energy than the active state, and if a metastable active state is present 
initially. These conditions lead to restrictions on the demand parameters a and h. If the demand 
is too large, the active state will immediately decay to the ground state and no metastable active 
state exists. If the demand is too small, the active state or some glassy state with lower activity 
is the ground state and no breakdowns can be observed. The region in which collapses can occur 
also depends on the topological exponent 7 of the network and on the delivery exponent /i. If 
is increased, the region in which avalanches can occur shrinks gradually. 

Random malfunction in the suppliers is modeled by a temperature parameter T. At finite 
temperatures, thermal fluctuations can cause network collapses, which arc prevented at T = by 
the metastability of the active state. During a collapse the roles of the hubs and of the poorly 
connected nodes has been monitored separately. It has been possible to identify which type of 
nodes is responsible for initiating the breakdown as a fimction of the distribution and network 
parameters. The lowest temperature T(, at which a network blackout takes place depends strongly 
on the different parameters in the model. Increasing this "breakdown temperature" is of great 
interest in the protection of the distribution system. The most stable situation appears to be that 
in which there are many large suppliers in the network. For large values of /x, this requires many 
poorly connected nodes while for /i = many hubs are needed. Rendering the amount of goods 
every consumer receives more homogeneous throughout the network also reduces the impact of 
thermal fluctuations in the system. Such a procedure can for instance be realized by decreasing 
Not all model parameters are independent. We have identified a scaled or effective temperature 
variable 6 which incorporates the 7 and fi dependencies in such a way that a good data collapse 
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can occur when the network activity for various cases is measured as a function of 6 instead of 
T. While this scahng is very useful at the level of a mean-field approximation, we also observed 
that fluctuations in the node degree, which are large for networks with 7 < 3, are responsible for 
deviations from simple scaling if also 7' = (7 — /Lt)/(1 — /i) < 3. 

The similarity between the model and real-life distribution systems can be improved in various 
ways. In the distribution model all the consumers depend on two suppliers. Using a bipartite 
network, i.e., a network with two kinds of nodes with links running only between unlike nodes, 
we can extend the model to incorporate an unrestricted number of suppliers for each consumer, 
reflecting the consumer's freedom of choice. Another extension is concerned with partly active 
suppliers. In the current model, a node is active or inactive. Using continuous spins, or discrete 
Potts spins, we could also model suppliers with tunable activity. Extensions of our model could 
also describe systems that evolve in time, for instance by implementing the model on growing 
networks or on networks in which rewiring is possible. Time-dependent demand parameters a and 
h could be used to model the evolving economic characteristics of the consumers, etc. 

We conclude that the distribution model introduced in this Paper offers a possible starting 
point for studying collapses in certain real-life distribution systems, from the viewpoint of the 
phenomenology of dynamical critical phenomena with a non-conserved order parameter. Note 
that in contrast with (most) other models of distribution or transportation networks [22] there is 
no conservation law or continuity equation in our network. The amounts of goods flowing along 
the edges of the network arc stochastic variables controlled by fluctuating spin states. In this sense 
our distribution network based on spin variables provides a complementary approach. 
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Figure 1: Phase diagram in (a, &)-space of the distribution model for a scale-free network with 1000 
nodes, topological exponent 7 = 3 and minimal degree m = 2. Part (a) corresponds to = 0, 
while part (b) pertains to /i = 0.2. The uppermost (dotted) lines indicate the upper boundaries 
of metastable active states for single-spin-flip dynamics at T = 0. The second uppermost lines 
(solid) denote the absolute stability limit of the inactive state at T = 0. The dark-grey region 
filling the space between the uppermost and second uppermost lines corresponds to a regime in 
which, at T > 0, a single-spin flip can cause an avalanche starting from a metastable active state. 
The next lower line (dash-dotted; coincident with the foregoing one in part (a) but not in part 
(b) of the figure) marks the absolute stability limit of the active state at T = 0. In Fig. (b), the 
ground state is unknown in the no-man's land between the two absolute stability lines, indicated by 
the light-grey shading. Finally the lowermost lines delineates the lower boundaries of metastable 
inactive states for single-spin- flip dynamics at T = 0. In part (b) the white cross indicates the 
specific values of a and b used in Fig. 3(a). 
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Figure 2: (a) Comparison of two mean- field and two exact phase boundaries for a distribution 
model with topological exponent 7 = 3 and delivery exponent = 0.2, with 1000 nodes and 
m = 2. Shown are the upper limit of the metastable active state for single-spin- flip dynamics 
at T = (uppermost curves; dotted line for the exact result and dashed line for the mean-field 
approximation) and the upper limit of absolute stability of the active state at T = (lowermost 
curves; dash-dotted line for the exact result and dashed line for the mean-field approximation); 
(b) Phase diagram for 7 = 3 and ^ — 1. The uppermost broken curve is the upper limit of 
metastability of the active state for single-spin flips at T = in the mean-field approximation 
(dashed line). The smooth curve below it (dotted) is the numerically exact result for that same 
metastability limit. Also shown are two crossing lines. The upper segments (solid) form the lower 
limit of absolute stability of the inactive state at T = 0, while the lower segments (dash-dotted) 
mark the upper limit of absolute stability of the active state at T = 0. The region that is (most 
likely) susceptible of avalanches is the small triangle shaded in dark grey. The no-man's land, with 
an unknown ground state, is shaded in light grey. The white crosses indicate the specific values of 
a and b used in Figs. 3(b-d). 
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Figure 3: (Color online) Typical examples of breakdowns in distribution networks. The graphs 
show the mean magnetization M in the network as a function of time. All figures are for the 
distribution model on scale-free networks with 7 = 3 and 1000 nodes. The values of the constants 
a, b and differ in the different subfigures. They correspond to the positions of the crosses in the 
phase diagrams Fig. 1(b) and Fig. 2(b) for /i = 0.2 and /x = 1, respectively. The solid (black; 
middle) curve indicates the mean magnetization, averaged over all possible node degrees. The dot- 
dashed (green; top in (a) and (b), bottom in (c) and (d)) curve gives the mean magnetization of 
the nodes with 8 or more neighbours {M{q > 8)) while the dashed (red; bottom in (a) and (b), top 
in (c) and (d)) curve shows the mean magnetization of the nodes with 2 neighbours {M{q — 2)). 
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Figure 4: (Color online) Effect of the delivery exponent fi on the thermal fluctuations. The 
figures show the "final" mean magnetization, obtained after 4000 time steps, as a function of the 
temperature T (Figs. 4(a) and 4(c)) and as a function of the effective or scaled temperature O 
(Figs. 4(b) and 4(d)) for different values of /i. All simulations are done on scale-free networks with 
1000 nodes, m = 2 and 7 = 3. We used a = 1.0 and b = 0.4 in Figs. 4(a) and 4(b), while a = 0.85 
and b = 0.4 in Figs. 4(c) and 4(d). All data points are averages over 10 network realizations. 
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Figure 5: (Color online) Effect of the topology on the thermal fluctuations. The figures show 
the magnetization as a function of the temperature T (Figs, (a) and (c)) and the parameter Q 
(Figs, (b) and (d)) for different values of 7. We study two different regimes depending on the 
parameter fj,. In Figs, (a) and (b), fi — 0, while /i = 1 in Figs, (c) and (d). All simulations were 
done on scale-free networks with 1000 nodes, with demand constants a = 1.0 and b = 0.2. All data 
points are averages over 10 network realizations. 



22 



